function c = chebft(a, b, n, tau)
% chebyshev fit; from Press et al, p 192.
% matlab code P. Manis 21 Mar 1999
% "Given a function func, lower and upper limits of the interval [a,b], and
% a maximum degree, n, this routine computes the n coefficients c[1..n] such that
% func(x) sum(k=1, n) of ck*Tk(y) - c0/2, where y = (x -0.5*(b+a))/(0.5*(b-a))
% This routine is to be used with moderately large n (30-50) the array of c's is 
% subsequently truncated at the smaller value m such that cm and subsequent
% terms are negligible."
%
bma = 0.5*(b-a);
bpa = 0.5*(b+a);
for k=0:n-1
   y = cos(pi*(k+0.5)/n);
 %  disp(sprintf('y=%8.3f, x=%8.3f', y, y*bma+bpa))
   f(k+1) = 1*exp(-(y*bma+bpa)/tau);
%   disp(sprintf('x@ = %8.3f', y*bma+bpa))
end
%fac = 2/n;
fac = 1;
for j=0:n-1
   sum=0;
   sum2 = 0;
   for k=0:n-1
      R = cos(pi*j*(k+0.5)/n);
      sum = sum+ f(k+1)*R;
      sum2 = sum2 + R^2;
   end
   c(j+1)=fac*sum/sum2;
end


